Analysis only KO Control

library(edgeR)
library (lattice)
library (snpStats)
library (tidyverse)
library(EnhancedVolcano)
library(VennDiagram)
library(gplots)
library(dplyr)

library(fgsea)
library(data.table)
library(ggplot2)
library(Rcpp)

counts = read.delim('counts_matrix_KO_Control.csv', row.names = 1, sep = ',')
head(counts)

names(counts)

[1] "D0.12.3" "D0.8.1"  "D0.28.3" "D0.12.1" "D0.G.1"  "D0.G.3"  "D0.G.2" 

d0 = DGEList(counts)
d0 = calcNormFactors(d0)
head(d0)

An object of class "DGEList"
$counts
                  D0.12.3 D0.8.1 D0.28.3 D0.12.1 D0.G.1 D0.G.3 D0.G.2
ENSG00000223972.5       2      3       2       4      4      4      3
ENSG00000227232.5     200    136     151     195    199    174    195
ENSG00000278267.1       0      0       0       0      0      0      0
ENSG00000243485.5       0      0       0       0      0      0      0
ENSG00000284332.1       0      0       0       0      0      0      0
ENSG00000237613.2       0      0       0       0      0      0      0

$samples

cutoff = 1
drop = which(apply(cpm(d0), 1, max) < cutoff)
d = d0[-drop,] 
dim(d) # number of genes left

[1] 14354     7

print(d)

An object of class "DGEList"
$counts
                  D0.12.3 D0.8.1 D0.28.3 D0.12.1 D0.G.1 D0.G.3 D0.G.2
ENSG00000227232.5     200    136     151     195    199    174    195
ENSG00000279457.4     284    282     208     225    284    312    228
ENSG00000225972.1    1128    643     802    1167    784    734    965
ENSG00000225630.1    1928   1398    1521    2028   1442   1637   1859
ENSG00000237973.1   13127   9470   10499   14625  10132  11545  12922
14349 more rows ...

$samples

NA

sample_description = read.delim('sample_descriptions_KO_Control.csv', sep = ',')
sample_description

#denoted tb, group treatment and background. fulv.Control or palbo.KO
group_treatment_background=interaction(sample_description$Treatment,sample_description$Background)
print(group_treatment_background)

[1] Pretx.KO      Pretx.KO      Pretx.KO      Pretx.KO      Pretx.Control Pretx.Control Pretx.Control
Levels: Pretx.Control Pretx.KO

plotMDS(d, col = as.numeric(group_treatment_background))

#denoted tdb, group treatment and days and background. fulv.0.control or palbo.16.KO
group_treatment_days_background = interaction(sample_description$Treatment,sample_description$Days,sample_description$Background)
print(group_treatment_days_background)

[1] Pretx.0.KO      Pretx.0.KO      Pretx.0.KO      Pretx.0.KO      Pretx.0.Control Pretx.0.Control Pretx.0.Control
Levels: Pretx.0.Control Pretx.0.KO

plotMDS(d, col = as.numeric(group_treatment_days_background))

#denoted tdg, group treatment and days and guide. fulv.0.12 or palbo.16.G
group_treatment_days_guide = interaction(sample_description$Treatment,sample_description$Days,sample_description$KO_Guide)
print(group_treatment_days_guide)

[1] Pretx.0.12 Pretx.0.8  Pretx.0.28 Pretx.0.12 Pretx.0.G  Pretx.0.G  Pretx.0.G 
Levels: Pretx.0.12 Pretx.0.28 Pretx.0.8 Pretx.0.G

plotMDS(d, col = as.numeric(group_treatment_days_guide))

NA
NA
NA
NA

grouping_tdb = group_treatment_days_background


mm_tdb = model.matrix(~0 + grouping_tdb)
y_tdb = voom(d, mm_tdb, plot = T)

fit_tdb = lmFit(y_tdb,mm_tdb)
head(fit_tdb)

An object of class "MArrayLM"
$coefficients
                  grouping_tdbPretx.0.Control grouping_tdbPretx.0.KO
ENSG00000227232.5                    1.756831               1.581529
ENSG00000279457.4                    2.284249               2.155806
ENSG00000225972.1                    3.874888               4.002993
ENSG00000225630.1                    4.869605               4.916537
ENSG00000237973.1                    7.678326               7.719442
ENSG00000229344.1                    6.304919               6.347772

$stdev.unscaled
                  grouping_tdbPretx.0.Control grouping_tdbPretx.0.KO
ENSG00000227232.5                  0.12280460             0.11190515
ENSG00000279457.4                  0.10390522             0.09382165
ENSG00000225972.1                  0.06337341             0.05251130
ENSG00000225630.1                  0.04837761             0.04134976
ENSG00000237973.1                  0.03181627             0.02748405
ENSG00000229344.1                  0.03663558             0.03156744

$sigma
[1] 1.370034 1.173333 4.041833 3.805404 6.225918 4.988128

$df.residual
[1] 5 5 5 5 5 5

$cov.coefficients
                            grouping_tdbPretx.0.Control grouping_tdbPretx.0.KO
grouping_tdbPretx.0.Control                   0.3333333                   0.00
grouping_tdbPretx.0.KO                        0.0000000                   0.25

$pivot
[1] 1 2

$rank
[1] 2

$Amean
ENSG00000227232.5 ENSG00000279457.4 ENSG00000225972.1 ENSG00000225630.1 ENSG00000237973.1 ENSG00000229344.1 
         1.669985          2.212636          3.965151          4.907124          7.705747          6.335550 

$method
[1] "ls"

$design
  grouping_tdbPretx.0.Control grouping_tdbPretx.0.KO
1                           0                      1
2                           0                      1
3                           0                      1
4                           0                      1
5                           1                      0
6                           1                      0
7                           1                      0
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$grouping_tdb
[1] "contr.treatment"

print(colnames(coef(fit_tdb)))

[1] "grouping_tdbPretx.0.Control" "grouping_tdbPretx.0.KO"     

contrast_group_ko_c= makeContrasts(grouping_tdbPretx.0.KO - grouping_tdbPretx.0.Control, levels = colnames(coef(fit_tdb)))

#contrast
contrastFunction = function(fit_object, contrast_object){
  tmp = contrasts.fit(fit_object, contrast_object)
  tmp = eBayes(tmp)
  top.table = topTable(tmp, sort.by = 'logFC', n = Inf)
  head(top.table, 50)
  qqunif.plot(top.table$P.Val)
  hist(top.table$P.Val)
  print(paste0('number gene adj.p.val <.05 = ', length(which(top.table$adj.P.Val <.05))))

  df = as.data.frame(top.table)
  
  #replace add genename to esng
  gene_name_df <- read.csv('name_description.csv', sep = ',')

  df$gene_name = gene_name_df$Description[match(row.names(df), gene_name_df$Name)]
  return(df)
}



#qqplot
qqunif.plot<-function(pvalues, 
    should.thin=T, thin.obs.places=2, thin.exp.places=2, 
    xlab=expression(paste("Expected (",-log[10], " p-value)")),
    ylab=expression(paste("Observed (",-log[10], " p-value)")), 
    draw.conf=TRUE, conf.points=1000, conf.col="lightgray", conf.alpha=.05,
    already.transformed=FALSE, pch=20, aspect="iso", prepanel=prepanel.qqunif,
    par.settings=list(superpose.symbol=list(pch=pch)), ...) {
    
    
    #error checking
    if (length(pvalues)==0) stop("pvalue vector is empty, can't draw plot")
    if(!(class(pvalues)=="numeric" || 
        (class(pvalues)=="list" && all(sapply(pvalues, class)=="numeric"))))
        stop("pvalue vector is not numeric, can't draw plot")
    if (any(is.na(unlist(pvalues)))) stop("pvalue vector contains NA values, can't draw plot")
    if (already.transformed==FALSE) {
        if (any(unlist(pvalues)==0)) stop("pvalue vector contains zeros, can't draw plot")
    } else {
        if (any(unlist(pvalues)<0)) stop("-log10 pvalue vector contains negative values, can't draw plot")
    }
    
    
    grp<-NULL
    n<-1
    exp.x<-c()
    if(is.list(pvalues)) {
        nn<-sapply(pvalues, length)
        rs<-cumsum(nn)
        re<-rs-nn+1
        n<-min(nn)
        if (!is.null(names(pvalues))) {
            grp=factor(rep(names(pvalues), nn), levels=names(pvalues))
            names(pvalues)<-NULL
        } else {
            grp=factor(rep(1:length(pvalues), nn))
        }
        pvo<-pvalues
        pvalues<-numeric(sum(nn))
        exp.x<-numeric(sum(nn))
        for(i in 1:length(pvo)) {
            if (!already.transformed) {
                pvalues[rs[i]:re[i]] <- -log10(pvo[[i]])
                exp.x[rs[i]:re[i]] <- -log10((rank(pvo[[i]], ties.method="first")-.5)/nn[i])
            } else {
                pvalues[rs[i]:re[i]] <- pvo[[i]]
                exp.x[rs[i]:re[i]] <- -log10((nn[i]+1-rank(pvo[[i]], ties.method="first")-.5)/(nn[i]+1))
            }
        }
    } else {
        n <- length(pvalues)+1
        if (!already.transformed) {
            exp.x <- -log10((rank(pvalues, ties.method="first")-.5)/n)
            pvalues <- -log10(pvalues)
        } else {
            exp.x <- -log10((n-rank(pvalues, ties.method="first")-.5)/n)
        }
    }


    #this is a helper function to draw the confidence interval
    panel.qqconf<-function(n, conf.points=1000, conf.col="gray", conf.alpha=.05, ...) {
        require(grid)
        conf.points = min(conf.points, n-1);
        mpts<-matrix(nrow=conf.points*2, ncol=2)
            for(i in seq(from=1, to=conf.points)) {
                    mpts[i,1]<- -log10((i-.5)/n)
                    mpts[i,2]<- -log10(qbeta(1-conf.alpha/2, i, n-i))
                    mpts[conf.points*2+1-i,1]<- -log10((i-.5)/n)
                    mpts[conf.points*2+1-i,2]<- -log10(qbeta(conf.alpha/2, i, n-i))
            }
            grid.polygon(x=mpts[,1],y=mpts[,2], gp=gpar(fill=conf.col, lty=0), default.units="native")
        }

    #reduce number of points to plot
    if (should.thin==T) {
        if (!is.null(grp)) {
            thin <- unique(data.frame(pvalues = round(pvalues, thin.obs.places),
                exp.x = round(exp.x, thin.exp.places),
                grp=grp))
            grp = thin$grp
        } else {
            thin <- unique(data.frame(pvalues = round(pvalues, thin.obs.places),
                exp.x = round(exp.x, thin.exp.places)))
        }
        pvalues <- thin$pvalues
        exp.x <- thin$exp.x
    }
    gc()
    
    prepanel.qqunif= function(x,y,...) {
        A = list()
        A$xlim = range(x)*1.02
        A$xlim[1]=0
        A$ylim = range(y)*1.02
        A$ylim[1] = 0
        return(A)
    }

    #draw the plot
    xyplot(pvalues~exp.x, groups=grp, xlab=xlab, ylab=ylab, aspect=aspect,
        prepanel=prepanel, scales=list(axs="i"), pch=pch,
        panel = function(x, y, ...) {
            if (draw.conf) {
                panel.qqconf(n, conf.points=conf.points, 
                    conf.col=conf.col, conf.alpha=conf.alpha)
            };
            panel.xyplot(x,y, ...);
            panel.abline(0,1);
        }, par.settings=par.settings, ...
    )
}

#volcano plot
plotVolcano = function(df, select_lab, title= NULL,keyvals= NULL){
  EnhancedVolcano(df, lab = df$gene_name, selectLab = select_lab, x = 'logFC', y = 'adj.P.Val',xlim = c(min(df[['logFC']], na.rm = TRUE) - .5, max(df[['logFC']], na.rm = TRUE) +
    .5),
  ylim = c(0, max(-log10(df[['adj.P.Val']]), na.rm = TRUE) + .5),FCcutoff = .001,pCutoff = 1, cutoffLineType = 'blank' ,legendPosition = 'right',drawConnectors =TRUE,widthConnectors = 0.3, pointSize = 1.5, labSize = 3.0, maxoverlapsConnectors = 50, title = title,colCustom = keyvals )
}


gsea_res = function(ranked_df_object, pathway_object){
  set.seed(42)
  fgseaRes = fgsea(pathways = pathway_object, stats = ranked_df_object, minSize = 15, maxSize = 500, eps = 0)
  topPathwaysUp = fgseaRes[ES > 0][head(order(pval), n=Inf),]
  topPathwaysDown = fgseaRes[ES < 0][head(order(pval), n=Inf),]
  topPathwaysUp_pathways = fgseaRes[ES > 0][head(order(pval), n=10),pathway]
  topPathwaysDown_pathways = fgseaRes[ES < 0][head(order(pval), n=10),pathway]
  topPathwaysAll = c(topPathwaysUp_pathways, rev(topPathwaysDown_pathways))
  
  new_list = list(topPathwaysUp, topPathwaysDown, topPathwaysAll, fgseaRes)
  return(new_list)
    
}

df_tdb_ko_c = contrastFunction(fit_tdb, contrast_group_ko_c)

[1] "number gene adj.p.val <.05 = 1978"

#remove duplicate gene values (most genes have similar logFC. remove lower logFC. for gene CYB561D2, logFC is .2 and -.0711, for RGS5, logFC is -.1461 and .004213093)
#duplicates= df_tdb_ko_c[duplicated(df_tdb_ko_c$gene_name),]
df_tdb_ko_c = df_tdb_ko_c[!duplicated(df_tdb_ko_c$gene_name),]

#save as csv
#write.csv(df_tdb_ko_c, 'df_tdb_ko_c.csv')

plotVolcano(df = df_tdb_ko_c, select_lab = df_tdb_ko_c$gene_name)

Warning: ggrepel: 14297 unlabeled data points (too many overlaps). Consider increasing max.overlaps

#t statistic for gsea
ranked_df = df_tdb_ko_c
ranked_df = within(ranked_df, rank <- ranked_df$t)
ranked_df = ranked_df[,c('gene_name', 'rank')]
ranked_df = ranked_df %>% arrange(rank)
ranked_df = deframe(ranked_df)




#plotGseaTable(h.all_pathway[topPathways_h_all_df], ranked_df, fgseaRes_h_all_df, gseaParam= 1, colwidths = c(17, 3, 2, 2, 2.5))

#c1 hallmark gene set
h.all_pathway = gmtPathways('h.all.v7.4.symbols.gmt.txt')
results_list = gsea_res(ranked_df, pathway_object = h.all_pathway)

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
  There are ties in the preranked stats (0.01% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

h.all_pathwayUp = results_list[1]
h.all_pathwayDown = results_list[2]
h.all_pathwayAll = results_list[3]
h.all_pathway_fgseaRes = results_list[4]
print(h.all_pathwayUp)

[[1]]

print(h.all_pathwayDown)

[[1]]

#plotGseaTable(h.all_pathway[h.all_pathwayAll, ranked_df, h.all_pathway_fgseaRes, gseaParam= 1, colwidths = c(17, 3, 2, 2, 2.5))

#c2 curated gene sets
c2.all_pathway = gmtPathways('c2.all.v7.4.symbols.gmt.txt')
results_list = gsea_res(ranked_df, pathway_object = c2.all_pathway)

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
  There are ties in the preranked stats (0.01% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

c2.all_pathwayUp = results_list[1]
c2.all_pathwayDown = results_list[2]
c2.all_pathwayAll = results_list[3]
c2.all_pathway_fgseaRes = results_list[4]
print(c2.all_pathwayUp)

[[1]]

print(c2.all_pathwayDown)

[[1]]

NA

#c3 regulatory target transcription factor targets
c3.tft_pathway = gmtPathways('c3.tft.v7.4.symbols.gmt.txt')
results_list = gsea_res(ranked_df, pathway_object = c3.tft_pathway)

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
  There are ties in the preranked stats (0.01% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

c3.tft_pathwayUp = results_list[1]
c3.tft_pathwayDown = results_list[2]
c3.tft_pathwayAll = results_list[3]
c3.tft_pathway_fgseaRes = results_list[4]
print(c3.tft_pathwayUp)

[[1]]

print(c3.tft_pathwayDown)

[[1]]

NA

#c5 gene ontology molecular function
c5.go.mf_pathway = gmtPathways('c5.go.mf.v7.4.symbols.gmt.txt')
results_list = gsea_res(ranked_df, pathway_object = c5.go.mf_pathway)

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
  There are ties in the preranked stats (0.01% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

c5.go.mf_pathwayUp = results_list[1]
c5.go.mf_pathwayDown = results_list[2]
c5.go.mf_pathwayAll = results_list[3]
c5.go.mf_pathway_fgseaRes = results_list[4]
print(c5.go.mf_pathwayUp)

[[1]]

print(c5.go.mf_pathwayDown)

[[1]]

NA

#c6 ongogenic signature genesets
c6.all_pathway = gmtPathways('c6.all.v7.4.symbols.gmt.txt')
results_list = gsea_res(ranked_df, pathway_object = c6.all_pathway)

Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
  There are ties in the preranked stats (0.01% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.

c6.all_pathwayUp = results_list[1]
c6.all_pathwayDown = results_list[2]
c6.all_pathwayAll = results_list[3]
c6.all_pathway_fgseaRes = results_list[4]
print(c6.all_pathwayUp)

[[1]]

print(c6.all_pathwayDown)

[[1]]

NA

#highlight specific genes in volcano plot
#gene_df<- read.csv("geneset_reg_esr1.txt")
df = df_tdb_ko_c
gene_df = h.all_pathway
df_gene_interest = df[(df$gene_name %in% gene_df$HALLMARK_ESTROGEN_RESPONSE_EARLY),]

#make new column with only labels for values in gene_list. use if gene_lsit is list
df$gene_label_interest = df_gene_interest$gene_name[match(df$gene_name, df_gene_interest$gene_name)]

keyvals = ifelse(is.na(df$gene_label_interest), 'royalblue','red')



names(keyvals)[keyvals == 'red'] <- 'interest'
names(keyvals)[keyvals == 'royalblue'] <- 'notinterest'
  
pointsize = c(ifelse(is.na(df$gene_label_interest),.5, 3))
plotVolcano(df = df, select_lab = df$gene_label_interest, title = 'Hallmark_estrogen_response_early',keyvals = keyvals)

Warning: ggrepel: 143 unlabeled data points (too many overlaps). Consider increasing max.overlaps

EnhancedVolcano(df_gene_interest, lab = df_gene_interest$gene_name, x = 'logFC', y = 'adj.P.Val',FCcutoff = .001,pCutoff = 1, cutoffLineType = 'blank' ,legendPosition = 'right',drawConnectors =FALSE,widthConnectors = 0.3, pointSize = 1.5, labSize = 3.0, maxoverlapsConnectors = 50)

Warning: Ignoring unknown parameters: xlim, ylim

plotEnrichment(h.all_pathway[['HALLMARK_ESTROGEN_RESPONSE_EARLY']],ranked_df) + labs(title = 'HALLMARK_ESTROGEN_RESPONSE_EARLY')

enrichmentplot = function(gsea_pathway, specific_pathway, ranked_dataframe = ranked_df){
  plotEnrichment(gsea_pathway[[specific_pathway]], ranked_dataframe) + labs(title = specific_pathway)
}

# interest_pathway = h.all_pathwayDown[[1]]$pathway
# enrichmentplot(gsea_pathway = h.all_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df)

for (i in 1:10){

  interest_pathway = h.all_pathwayUp[[1]]$pathway[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = h.all_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "HALLMARK_INTERFERON_GAMMA_RESPONSE"
[1] "HALLMARK_INTERFERON_ALPHA_RESPONSE"
[1] "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"
[1] "HALLMARK_ALLOGRAFT_REJECTION"
[1] "HALLMARK_MYOGENESIS"
[1] "HALLMARK_XENOBIOTIC_METABOLISM"
[1] "HALLMARK_UV_RESPONSE_UP"
[1] "HALLMARK_OXIDATIVE_PHOSPHORYLATION"
[1] "HALLMARK_COMPLEMENT"
[1] "HALLMARK_COAGULATION"

for (i in 1:5){

  interest_pathway = h.all_pathwayDown[[1]]$pathway[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = h.all_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "HALLMARK_ESTROGEN_RESPONSE_EARLY"
[1] "HALLMARK_ESTROGEN_RESPONSE_LATE"
[1] "HALLMARK_G2M_CHECKPOINT"
[1] "HALLMARK_MITOTIC_SPINDLE"
[1] "HALLMARK_ANDROGEN_RESPONSE"

c2_up_down_interest = read.csv('c2_up_down_pathways.csv')


c2_interest_list_Up = c2_up_down_interest$c2.up

for (i in 1:length(c2_interest_list_Up)){

  interest_pathway = c2_interest_list_Up[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = c2.all_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_DN"
[1] "CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_DN"
[1] "LIM_MAMMARY_STEM_CELL_UP"
[1] "BERTUCCI_MEDULLARY_VS_DUCTAL_BREAST_CANCER_UP"
[1] "SCHUETZ_BREAST_CANCER_DUCTAL_INVASIVE_UP"
[1] "STEIN_ESRRA_TARGETS_UP"
[1] "SMID_BREAST_CANCER_NORMAL_LIKE_UP"
[1] "SMID_BREAST_CANCER_LUMINAL_B_DN"
[1] "CHICAS_RB1_TARGETS_CONFLUENT"
[1] "SANSOM_APC_TARGETS_DN"
[1] "WP_NUCLEAR_RECEPTORS_METAPATHWAY"
[1] "STEIN_ESRRA_TARGETS"
[1] "SATO_SILENCED_BY_METHYLATION_IN_PANCREATIC_CANCER_1"
[1] "VERRECCHIA_EARLY_RESPONSE_TO_TGFB1"
[1] "LIM_MAMMARY_LUMINAL_PROGENITOR_UP"
[1] "GOTZMANN_EPITHELIAL_TO_MESENCHYMAL_TRANSITION_UP"
[1] "HELLER_HDAC_TARGETS_SILENCED_BY_METHYLATION_UP"
[1] "HUANG_FOXA2_TARGETS_DN"
[1] "SMID_BREAST_CANCER_BASAL_UP"
[1] "JECHLINGER_EPITHELIAL_TO_MESENCHYMAL_TRANSITION_UP"
[1] "BOWIE_RESPONSE_TO_TAMOXIFEN"
[1] "GOZGIT_ESR1_TARGETS_UP"
[1] "VANTVEER_BREAST_CANCER_ESR1_DN"

c2_interest_list_Down = c2_up_down_interest$c2.down

for (i in 1:length(c2_interest_list_Down)){

  interest_pathway = c2_interest_list_Down[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = c2.all_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "CHARAFE_BREAST_CANCER_LUMINAL_VS_MESENCHYMAL_UP"
[1] "CREIGHTON_ENDOCRINE_THERAPY_RESISTANCE_5"
[1] "CHARAFE_BREAST_CANCER_LUMINAL_VS_BASAL_UP"
[1] "FARMER_BREAST_CANCER_BASAL_VS_LULMINAL"
[1] "MASSARWEH_RESPONSE_TO_ESTRADIOL"
[1] "VANTVEER_BREAST_CANCER_ESR1_UP"
[1] "PUJANA_BRCA_CENTERED_NETWORK"
[1] "MASSARWEH_TAMOXIFEN_RESISTANCE_DN"
[1] "REACTOME_CHROMATIN_MODIFYING_ENZYMES"
[1] "ZHANG_BREAST_CANCER_PROGENITORS_UP"
[1] "DUTERTRE_ESTRADIOL_RESPONSE_24HR_UP"
[1] "FISCHER_G1_S_CELL_CYCLE"
[1] "PUJANA_BRCA2_PCC_NETWORK"
[1] "DOANE_BREAST_CANCER_ESR1_UP"
[1] "MCBRYAN_PUBERTAL_BREAST_4_5WK_UP"
[1] "DUTERTRE_ESTRADIOL_RESPONSE_6HR_DN"
[1] "KONG_E2F3_TARGETS"
[1] "CREIGHTON_ENDOCRINE_THERAPY_RESISTANCE_1"
[1] "REACTOME_HATS_ACETYLATE_HISTONES"
[1] "HOLLERN_EMT_BREAST_TUMOR_DN"
[1] "CREIGHTON_ENDOCRINE_THERAPY_RESISTANCE_4"
[1] "EGUCHI_CELL_CYCLE_RB1_TARGETS"
[1] "REACTOME_HDMS_DEMETHYLATE_HISTONES"

for (i in 1:20){

  interest_pathway = c3.tft_pathwayUp[[1]]$pathway[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = c3.tft_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "NFE2_01"
[1] "GGGNNTTTCC_NFKB_Q6_01"
[1] "NRF2_Q4"
[1] "AP1_Q6_01"
[1] "SF1_Q6"
[1] "AP1_Q4_01"
[1] "AP1_C"
[1] "TGASTMAGC_NFE2_01"
[1] "ISRE_01"
[1] "SREBP1_Q6"
[1] "AP1_Q6"
[1] "E47_01"
[1] "NR0B1_TARGET_GENES"
[1] "TGACCTTG_SF1_Q6"
[1] "BACH2_01"
[1] "IRF_Q6"
[1] "STAT5B_01"
[1] "CCAWWNAAGG_SRF_Q4"
[1] "AP1_01"
[1] "RGAGGAARY_PU1_Q6"

for (i in 1:20){

  interest_pathway = c3.tft_pathwayDown[[1]]$pathway[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = c3.tft_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "HSD17B8_TARGET_GENES"
[1] "SMTTTTGT_UNKNOWN"
[1] "E2F_Q3"
[1] "E2F_Q3_01"
[1] "DMRT1_TARGET_GENES"
[1] "E2F1_Q4_01"
[1] "TGACATY_UNKNOWN"
[1] "TGTTTGY_HNF3_Q6"
[1] "AREB6_04"
[1] "TTANWNANTGGM_UNKNOWN"
[1] "FOXO3_01"
[1] "FOXO1_01"
[1] "E2F1_Q6_01"
[1] "AP3_Q6"
[1] "FAC1_01"
[1] "PITX2_Q2"
[1] "KCCGNSWTTT_UNKNOWN"
[1] "FOXO4_02"
[1] "SRY_01"
[1] "GGATTA_PITX2_Q2"

for (i in 1:20){

  interest_pathway = c5.go.mf_pathwayUp[[1]]$pathway[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = c5.go.mf_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "GOMF_ION_TRANSMEMBRANE_TRANSPORTER_ACTIVITY"
[1] "GOMF_ELECTRON_TRANSFER_ACTIVITY"
[1] "GOMF_PROTEIN_HOMODIMERIZATION_ACTIVITY"
[1] "GOMF_ACTIVE_TRANSMEMBRANE_TRANSPORTER_ACTIVITY"
[1] "GOMF_CALCIUM_ION_BINDING"
[1] "GOMF_MONOCARBOXYLIC_ACID_BINDING"
[1] "GOMF_ANION_TRANSMEMBRANE_TRANSPORTER_ACTIVITY"
[1] "GOMF_PEPTIDE_BINDING"
[1] "GOMF_ACTIVE_ION_TRANSMEMBRANE_TRANSPORTER_ACTIVITY"
[1] "GOMF_AMIDE_BINDING"
[1] "GOMF_CATION_TRANSMEMBRANE_TRANSPORTER_ACTIVITY"
[1] "GOMF_ABC_TYPE_TRANSPORTER_ACTIVITY"
[1] "GOMF_STRUCTURAL_MOLECULE_ACTIVITY"
[1] "GOMF_ORGANIC_ANION_TRANSMEMBRANE_TRANSPORTER_ACTIVITY"
[1] "GOMF_G_PROTEIN_COUPLED_RECEPTOR_BINDING"
[1] "GOMF_MOLECULAR_TRANSDUCER_ACTIVITY"
[1] "GOMF_PASSIVE_TRANSMEMBRANE_TRANSPORTER_ACTIVITY"
[1] "GOMF_SECONDARY_ACTIVE_TRANSMEMBRANE_TRANSPORTER_ACTIVITY"
[1] "GOMF_ENZYME_INHIBITOR_ACTIVITY"
[1] "GOMF_OXIDOREDUCTASE_ACTIVITY_ACTING_ON_NAD_P_H"

for (i in 1:20){

  interest_pathway = c5.go.mf_pathwayDown[[1]]$pathway[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = c5.go.mf_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "GOMF_DNA_BINDING_TRANSCRIPTION_ACTIVATOR_ACTIVITY"
[1] "GOMF_CHROMATIN_BINDING"
[1] "GOMF_DNA_BINDING_TRANSCRIPTION_REPRESSOR_ACTIVITY"
[1] "GOMF_CHROMATIN_DNA_BINDING"
[1] "GOMF_PROTEIN_DEMETHYLASE_ACTIVITY"
[1] "GOMF_MRNA_BINDING"
[1] "GOMF_DNA_BINDING_TRANSCRIPTION_FACTOR_BINDING"
[1] "GOMF_ACTIN_BINDING"
[1] "GOMF_NUCLEOSOME_BINDING"
[1] "GOMF_DNA_SECONDARY_STRUCTURE_BINDING"
[1] "GOMF_LYS48_SPECIFIC_DEUBIQUITINASE_ACTIVITY"
[1] "GOMF_FOUR_WAY_JUNCTION_DNA_BINDING"
[1] "GOMF_CYTOKINE_RECEPTOR_ACTIVITY"
[1] "GOMF_RNA_POLYMERASE_CORE_ENZYME_BINDING"
[1] "GOMF_DEMETHYLASE_ACTIVITY"
[1] "GOMF_TRANSCRIPTION_COREGULATOR_ACTIVITY"
[1] "GOMF_HISTONE_LYSINE_N_METHYLTRANSFERASE_ACTIVITY"
[1] "GOMF_RNA_POLYMERASE_II_SPECIFIC_DNA_BINDING_TRANSCRIPTION_FACTOR_BINDING"
[1] "GOMF_IMMUNE_RECEPTOR_ACTIVITY"
[1] "GOMF_PEPTIDE_HORMONE_BINDING"

for (i in 1:20){

  interest_pathway = c6.all_pathwayUp[[1]]$pathway[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = c6.all_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "LTE2_UP.V1_DN"
[1] "RAF_UP.V1_UP"
[1] "KRAS.DF.V1_UP"
[1] "VEGF_A_UP.V1_UP"
[1] "EGFR_UP.V1_UP"
[1] "PGF_UP.V1_DN"
[1] "LEF1_UP.V1_UP"
[1] "NFE2L2.V2"
[1] "PTEN_DN.V1_DN"
[1] "PRC2_EZH2_UP.V1_UP"
[1] "KRAS.600.LUNG.BREAST_UP.V1_UP"
[1] "KRAS.LUNG.BREAST_UP.V1_UP"
[1] "BMI1_DN.V1_UP"
[1] "MEK_UP.V1_UP"
[1] "MEL18_DN.V1_UP"
[1] "JNK_DN.V1_DN"
[1] "ATF2_S_UP.V1_DN"
[1] "KRAS.LUNG_UP.V1_UP"
[1] "MTOR_UP.N4.V1_UP"
[1] "KRAS.600_UP.V1_UP"

for (i in 1:20){

  interest_pathway = c6.all_pathwayDown[[1]]$pathway[i]
  print(interest_pathway)
  print(enrichmentplot(gsea_pathway = c6.all_pathway, specific_pathway = interest_pathway, ranked_dataframe = ranked_df))
}

[1] "RAF_UP.V1_DN"
[1] "P53_DN.V1_UP"
[1] "STK33_NOMO_DN"
[1] "ESC_V6.5_UP_LATE.V1_UP"
[1] "STK33_NOMO_UP"
[1] "PRC1_BMI_UP.V1_UP"
[1] "JAK2_DN.V1_DN"
[1] "PDGF_UP.V1_DN"
[1] "LEF1_UP.V1_DN"
[1] "EGFR_UP.V1_DN"
[1] "STK33_UP"
[1] "GCNP_SHH_UP_EARLY.V1_UP"
[1] "MTOR_UP.N4.V1_DN"
[1] "TBK1.DF_DN"
[1] "CAHOY_ASTROCYTIC"
[1] "BCAT.100_UP.V1_DN"
[1] "PTEN_DN.V2_DN"
[1] "ESC_V6.5_UP_EARLY.V1_UP"
[1] "TGFB_UP.V1_DN"
[1] "SINGH_KRAS_DEPENDENCY_SIGNATURE"